library(stress.test.plot.report)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.3.2

parameters

# tiltrisk_input_path <- here::here("workspace","tiltrisk","tiltrisk_alpha.csv")
tiltrisk_input_path <- here::here("workspace","tiltrisk","tiltrisk_nrisk.csv") 
#         Change_NPV_shock 
#         Change_NPV_with_ecosystem_cost
#         Change_NPV_with_ecosystem_social_cost

tiltrisk_trajectories_input_path <- here::here("workspace","tiltrisk","tiltrisk_trajectories.csv")

dataload

tiltrisk_df <- readr::read_csv(tiltrisk_input_path) %>%
  rename(activity_name=`Activity Name`)

## Rows: 211620 Columns: 35
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (17): companies_id, company_name, country, main_activity, clustered, act...
## dbl (18): Change_NPV_shock, Change_NPV_with_ecosystem_cost, Change_NPV_with_...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

tiltrisk_df <- tiltrisk_df %>%
      mutate(
        pd_difference = pd_shock - pd_baseline,
        crispy_perc_value_change=Change_NPV_with_physical_risk
      )

dataprep

get_tiltrisk_plot_df <- function(tiltrisk_df, group_cols){
    tiltrisk_plot_df <- tiltrisk_df %>% 
      group_by_at(c(group_cols, "term")) %>%
      summarise(
        crispy_perc_value_change = mean(crispy_perc_value_change),
        pd_baseline=mean(pd_baseline),
        pd_shock=mean(pd_shock),
        .groups="drop"
      ) %>%
      mutate(
        pd_difference = pd_shock - pd_baseline
      )
    return(tiltrisk_plot_df)
}

basic_group_cols <- c("run_id", "baseline_scenario", "shock_scenario", "country")

tiltrisk_plot_df_companies_main_activity <- get_tiltrisk_plot_df(tiltrisk_df, c(basic_group_cols, "main_activity"))

tiltrisk_plot_df_companies_clustered <- get_tiltrisk_plot_df(tiltrisk_df, c(basic_group_cols, "clustered"))


tiltrisk_cluster_activity <- get_tiltrisk_plot_df(tiltrisk_df, c("run_id", "baseline_scenario", "shock_scenario", "country", "main_activity", "clustered"))

tiltrisk_plot_df_companies_main_activity %>%
  readr::write_csv(here::here("workspace","tiltrisk","tiltrisk_agg_main_activity.csv"))

tiltrisk_plot_df_companies_main_activity %>%
  readr::write_csv(here::here("workspace","tiltrisk","tiltrisk_agg_clustered.csv"))

plots

for (run_id_i in unique(tiltrisk_plot_df_companies_main_activity$run_id)){
  plot_data <- tiltrisk_plot_df_companies_main_activity %>%
    filter(run_id==run_id_i)
  

  shock_scenario = unique(plot_data$shock_scenario)
  for (loc in unique(plot_data$country)){
    crispy_npv_change_plot <- pipeline_crispy_npv_change_plot(
      plot_data |> dplyr::filter(term==1, country==loc),
      x_var = "main_activity"
      ) + 
      ggplot2::ggtitle(paste(shock_scenario, '-',  loc))
    
      
      print(crispy_npv_change_plot)
    }
}

for (run_id_i in unique(tiltrisk_plot_df_companies_main_activity$run_id)){
  plot_data <- tiltrisk_plot_df_companies_main_activity %>%
    filter(run_id==run_id_i)
  

  shock_scenario = unique(plot_data$shock_scenario)
  for (loc in unique(plot_data$country)){
    
      pd_term_plot <- pipeline_crispy_pd_term_plot(
        crispy_data_agg = plot_data %>% filter(country==loc),
        facet_var = "main_activity"
      ) + 
        ggplot2::ggtitle(paste(shock_scenario, '-',  loc))
      
      print(pd_term_plot)
    }
}

Crispy sensitivity analysis

# Calculate the average ROC per company for the second plot
df_averages <- tiltrisk_df |>
  dplyr::group_by(
    company_name,
    country,
    run_id,
    main_activity,
    shock_scenario,
    shock_year
  ) |>
  dplyr::summarise(
    average_roc = mean(crispy_perc_value_change),
    average_pd_diff = mean(pd_difference)
  ) %>%
  rename(
    `Average PD Difference` = average_pd_diff,
    `Average NPV % change` = average_roc
  )

## `summarise()` has grouped output by 'company_name', 'country', 'run_id',
## 'main_activity', 'shock_scenario'. You can override using the `.groups`
## argument.

# Reshape data to long format
long_data <- df_averages %>%
  tidyr::pivot_longer(
    cols = c(`Average PD Difference`, `Average NPV % change`),
    names_to = "metric",
    values_to = "value"
  )

# Calculate means and confidence intervals
data_summary <- long_data %>%
  group_by(
    run_id,
    country,
    shock_scenario,
    shock_year,
    main_activity,
    metric
  ) %>%
  summarise(
    mean = mean(value),
    se = sd(value) / sqrt(n()),
    ci_upper = mean + qt(0.975, df = n() - 1) * se,
    ci_lower = mean - qt(0.975, df = n() - 1) * se
  ) %>%
  ungroup()

## Warning: There were 8 warnings in `summarise()`.
## The first warning was:
## ℹ In argument: `ci_upper = mean + qt(0.975, df = n() - 1) * se`.
## ℹ In group 1: `run_id = "32a7d72b-480d-4de3-8522-6bfb40cfe2a8"`, `country =
##   "austria"`, `shock_scenario = "IPR2023_FPS"`, `shock_year = 2030`,
##   `main_activity = "agent/ representative"`, `metric = "Average NPV % change"`.
## Caused by warning in `qt()`:
## ! NaNs produced
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 7 remaining warnings.

## `summarise()` has grouped output by 'run_id', 'country', 'shock_scenario',
## 'shock_year', 'main_activity'. You can override using the `.groups` argument.

for (run_id in unique(data_summary$run_id)) {
  plot_data <- data_summary %>% dplyr::filter(.data$run_id == .env$run_id)
  for (loc in unique(plot_data$country)){
  p1 <-
    ggplot(
      plot_data %>% filter(country==loc),
      aes(x = main_activity, y = mean, fill = metric)
    ) +
    geom_bar(stat = "identity", position = position_dodge()) +
    geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper),
      width = 0.2,
      position = position_dodge(0.9)
    ) +
    facet_wrap(~metric, scales = "free_y") +
    scale_fill_manual(values = c(
      "Average NPV % change" = "#5D9324",
      "Average PD Difference" = "#BAB6B5"
    )) +
    scale_y_continuous(labels = scales::label_percent(accuracy = 1)) +
    r2dii.plot::theme_2dii() +
    theme(
      axis.text.x = element_text(angle = 45, hjust = 1),
      legend.position = "none"
    ) +
    labs(y = "Average Value", x = "Business Unit") +
    ggtitle(paste(
      loc,
      "-",
      plot_data[1, "shock_scenario"] %>% pull()
      # , "- shock year :",
      # plot_data[1, "shock_year"] %>% pull()
    ))

  print(p1)
  }
}

PD distrib plots

for (run_id_i in unique(tiltrisk_plot_df_companies_main_activity$run_id)){
  data_plot <- tiltrisk_df %>%
    filter(run_id==run_id_i, term==1)
  
  filtering <- data_plot %>% 
    distinct(country ,main_activity ,companies_id) %>% 
    group_by(country, main_activity) %>% 
    summarise(nrow=n()) %>%
    filter(nrow>10) %>%
    select(-nrow) %>%
    ungroup()

  data_plot <- data_plot %>% inner_join(filtering)
  
  shock_scenario <- unique(data_plot$shock_scenario)
  
  density_plot <- make_density_plots(data_plot,
    numeric_values = "pd_shock",
    density_var = "country",
    group_variable="main_activity"
  ) +
    ggplot2::ggtitle(paste0("Distribution of PD at shock for ", shock_scenario, "shock scenario"))
  
  print(density_plot)
}

## `summarise()` has grouped output by 'country'. You can override using the
## `.groups` argument.
## Joining with `by = join_by(country, main_activity)`
## `summarise()` has grouped output by 'country'. You can override using the
## `.groups` argument.
## Joining with `by = join_by(country, main_activity)`

Heterogeneity in Transition Risk

for (run_id_i in unique(tiltrisk_df$run_id)) {
  data_plot <- tiltrisk_df |> 
    dplyr::filter(run_id == run_id_i)

  agg_analysis_data <- data_plot |>
    # dplyr::filter(.data$net_present_value_difference != 0) |>
    dplyr::select(.data$company_name, .data$crispy_perc_value_change, .data$pd_difference) |>
    dplyr::group_by(.data$company_name) |>
    dplyr::summarise(
      crispy_perc_value_change = mean(crispy_perc_value_change),
      pd_difference = mean(pd_difference),
      .groups = "drop"
    )

  # Sorting categories based on value1 in descending order
  plot_data <- agg_analysis_data |>
    # sample_frac(0.1) |>
    dplyr::arrange(dplyr::desc(.data$crispy_perc_value_change)) |>
    dplyr::mutate(company_name = factor(.data$company_name, levels = .data$company_name)) |>
    tidyr::pivot_longer(cols = c("crispy_perc_value_change", "pd_difference"), names_to = "variable", values_to = "value")


  # Plotting
  p1 <- ggplot(plot_data %>% filter(variable == "crispy_perc_value_change"), aes(x = factor(company_name), y = value, group = variable)) +
    geom_step(color = "#5D9324", size = 1) +
    scale_y_continuous(
      labels = scales::percent_format(accuracy = 1),
      breaks = scales::pretty_breaks(n = 5)
    ) +
    geom_hline(yintercept = 0, color = "lightgray", linetype = "dashed", size = 0.5) +
    r2dii.plot::theme_2dii() +
    theme(
      axis.text.x = element_blank(), # element_text(angle = 90, vjust = 0.5),
      axis.ticks.x = element_blank(),
      axis.title.y = element_text(size = 11),
      strip.background = element_blank(),
      strip.placement = "outside",
      legend.position = "none"
    ) +
    labs(x = NULL, y = NULL) +
    guides(fill = NULL) +
    ylab("Mean company percent value change")



  # Function to create bins every 10 observations
  bin_data <- function(data, bin_size) {
    data <- data %>%
      mutate(bin = (as.numeric(row_number()) - 1) %/% bin_size) %>%
      group_by(bin) %>%
      summarise(
        avg = mean(value),
        min = min(value),
        max = max(value)
      ) %>%
      ungroup()
    return(data)
  }

  # Bin data every 10 observations
  binned_data <- bin_data(plot_data %>% filter(variable == "pd_difference"), round(nrow(plot_data) / 100))

  # Create the plot
  p2 <- ggplot(binned_data, aes(x = factor(bin), y = avg)) +
    geom_col(fill = "#BAB6B5") +
    geom_errorbar(aes(ymin = min, ymax = max), width = 0.2) +
    scale_y_continuous(
      labels = scales::percent_format(accuracy = 1),
      breaks = scales::pretty_breaks(n = 5)
    ) +
    r2dii.plot::theme_2dii() +
    theme(
      axis.text.x = element_blank(),
      axis.ticks.x = element_blank(),
      axis.title.y = element_text(size = 11),
      strip.background = element_blank(),
      strip.placement = "outside"
    ) +
    labs(x = NULL, y = NULL) +
    guides(fill = NULL) +
    ylab("Mean climate Transition-related PD difference")

  le_plot <- cowplot::plot_grid(p1, p2, ncol = 1, align = "v")

  title <- cowplot::ggdraw() +
    cowplot::draw_label(
      paste(
        data_plot[1, "shock_scenario"] %>% pull(),
        " - shock year :",
        data_plot[1, "shock_year"] %>% pull()
      ),
      fontface = "bold",
      x = 0.5,
      hjust = 0.5
    ) +
    theme(
      # add margin on the left of the drawing canvas,
      # so title is aligned with left edge of first plot
      plot.margin = margin(0, 0, 0, 7)
    )

  le_plot <- cowplot::plot_grid(
    title, le_plot,
    ncol = 1,
    # rel_heights values control vertical title margins
    rel_heights = c(0.1, 1)
  )

  print(le_plot)
}

## Warning: Use of .data in tidyselect expressions was deprecated in tidyselect 1.2.0.
## ℹ Please use `"company_name"` instead of `.data$company_name`
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: Use of .data in tidyselect expressions was deprecated in tidyselect 1.2.0.
## ℹ Please use `"crispy_perc_value_change"` instead of
##   `.data$crispy_perc_value_change`
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: Use of .data in tidyselect expressions was deprecated in tidyselect 1.2.0.
## ℹ Please use `"pd_difference"` instead of `.data$pd_difference`
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Funding

EU LIFE Project Grant

Scientific Transition Risk Exercises for Stress tests & Scenario Analysis has received funding from the European Union’s Life programme under Grant No. LIFE21-GIC-DE-Stress under the LIFE-2021-SAP-CLIMA funding call.